eegmp_analysis/tools/[MEG]PLS/MEGPLS_PIPELINE_v2.02b/MEGPLS_COMPONENTS/MEGpipeline_MakePLSmat.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reads 4D-AFNI files to generate .mat files for PLS input.  %
% Inputs: 4D-AFNI Input, 3D-AFNI Mask, Seed & Behavior Data. %
% Last modified: Jan. 15, 2014                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2013-2014, Michael J. Cheung
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


function PLSmat = MEGpipeline_MakePLSmat(PLSmatSettings)

% Load PLSmat input & settings:
PLSmat     = PLSmatSettings;
name       = PLSmat.name;
time       = PLSmat.time;
paths      = PLSmat.paths;
Opt.format = 'vector';

% Clear existing Errorlog & Diary:
if exist('ErrorLog_MakePLSmat.txt', 'file')
    system('rm ErrorLog_MakePLSmat.txt');
end
if exist('Diary_MakePLSmat.txt', 'file')
    system('rm Diary_MakePLSmat.txt');
end

diary Diary_MakePLSmat.txt
ErrLog = fopen('ErrorLog_MakePLSmat.txt', 'a');



%========================================%
% CHECK FILES, DIMENSIONS & ORIENTATION: %
%========================================%

% Get dimensions & orientation from first .BRIK file:
[~, Brik, Info, ~] = BrikLoad(paths.Afni4D{1}{1,1}, Opt);

Brik		= double(Brik);
BrikDim4D   = size(Brik);
BrikDim3D   = BrikDim4D(1:3);
BrikOrient  = Info.ORIENT_SPECIFIC;
SampleInfo  = Info;


% Check AFNI files:
MissingFiles = 0;
MismatchDims = 0;

for g = 1:length(name.GroupID)
	for s = 1:length(name.SubjID{g})
		for c = 1:length(name.CondID)
			
            % Check if file exists (Since no LoadFTmat):
            if ~exist(paths.Afni4D{g}{s,c}, 'file')
                fprintf(ErrLog, 'ERROR: Input AFNI file does not exist: \n %s \n\n',...
                    paths.Afni4D{g}{s,c});
                
                MissingFiles = 1;
                continue;
            else
                [~, Brik, Info, ~] = BrikLoad(paths.Afni4D{g}{s,c}, Opt);
            end
			
            % Check dimensions:
			if ~isequal(size(Brik), BrikDim4D) || ~isequal(Info.ORIENT_SPECIFIC, BrikOrient)
                fprintf(ErrLog, ['ERROR: File has different dimensions or orientation:'...
                    '\n File: %s \n Dims: %s \n Orient: %s \n\n'],...
                    paths.Afni4D{g}{s,c}, num2str(size(Brik)), Info.Orientation(:,1));
                
				MismatchDims = 1;
			end
			
		end  % Cond
	end  % Subj
end  % Group

if MissingFiles == 1
    error('ERROR: Input AFNI .BRIK file(s) missing.')
end
if MismatchDims == 1
	error('ERROR: All AFNI files must have the same Dimensions & Orientation.')
end

% In case of 3D AFNI files:
if length(BrikDim4D) == 3
    BrikDim4D   = [BrikDim3D, 1];
    BrikNumTime = 1;
else
    BrikNumTime = BrikDim4D(4);
end

Brik = [];  % Free memory



%===================================%
% ACQUIRE WHOLE-BRAIN MASK INDICES: %
%===================================%

if strcmp(PLSmat.MaskType, 'BrainMask')
    
    
    %--- If generating a data-specific mask: ---%
    %-------------------------------------------%
    if strcmp(PLSmat.DataSpecificMask, 'yes')
        
        % Only keep "inside" voxels shared by certain % of the datasets.
        clear LoadData BrainMask NumDatasets
        
        BrainMask  = [];
        MaskErrors = 0;
        t = PLSmat.TimeIndices(1);  % Inside should be same across time anyways.
        
        for g = 1:length(name.GroupID)
            for s = 1:length(name.SubjID{g})
                for c = 1:length(name.CondID)
                    
                    LoadData = LoadFTmat(paths.NormSource{g}{s,c,t}, 'MakePLSmat');
                    
                    if isempty(LoadData)
                        MaskErrors = 1;
                        continue;
                    end
                    
                    if isempty(BrainMask)
                        BrainMask   = LoadData.inside;
                        NumDatasets = 1;
                        
                    else
                        if ~isequal(size(BrainMask), size(LoadData.inside))
                            fprintf(ErrLog, ['ERROR: For mask creation, this dataset does'...
                                'not have the same dimensions as the first dataset loaded.'...
                                '\n %s \n\n'], paths.NormSource{g}{s,c,t});
                            
                            MaskErrors = 1;
                            continue;
                        end
                        
                        BrainMask   = BrainMask + LoadData.inside;  % Add all inside data
                        NumDatasets = NumDatasets + 1;
                    end
                    
                    LoadData = [];  % Free memory
                
                end  % Cond
            end  % Subj
        end  % Group
        
        if MaskErrors == 0
            BrainMask = BrainMask ./ NumDatasets;  % Average to get %
            BrainMask(find(BrainMask < PLSmat.MaskKeepVoxThresh)) = 0;  % Set voxels under thresh to 0.
            BrainMask(find(BrainMask > 0)) = 1;
            BrainMask = logical(BrainMask);
            
            if ~isequal(size(BrainMask), BrikDim3D)
                fprintf(ErrLog, ['ERROR: Data-specific brain mask has different dimensions'...
                    'than the input 4D-AFNI files. \n\n']);
                error('ERROR: Data-specific brain mask has different dims than input 4D-AFNI files.');
            end
            
            BrainMaskFlat = reshape(BrainMask, 1, []);       % Flatten to 1D-array (applied to flat data).
            BrainMaskFlatIndices = find(BrainMaskFlat > 0);  % Linear voxel-indices of mask areas.
            
            BrainMaskFlat = [];  % Free memory
            
        else
            fprintf(ErrLog, ['ERROR: Failed to generate data-specific brain mask.'...
                '\n This may be due to mismatched dimensions in the NormSource files.'...
                '\n This may also be caused by missing or broken files.']);
            error('ERROR: Failed to generate data-specific brain mask. See ErrorLog.');
        end
        
    end

    
    %--- If selecting a brain-mask file: ---%
    %---------------------------------------%
    if strcmp(PLSmat.DataSpecificMask, 'no');
        
        % Check for brain mask:
        if ~exist(PLSmat.BrainMaskFile, 'file')
            fprintf(ErrLog, 'ERROR: Input Brain-Mask does not exist: \n %s \n\n',...
                PLSmat.BrainMaskFile);
            
            error('ERROR: Input Brain-Mask file could not be found.')
        end
        
        % Move selected brain mask to AnalysisID folder:
        [~, OrigMaskName, ~] = fileparts(PLSmat.BrainMaskFile);
        CopiedMaskFile = [paths.AnalysisID,'/MASKS_TEMPLATES/',OrigMaskName,'.BRIK'];
        
        if exist(CopiedMaskFile, 'file')
            delete(CopiedMaskFile);
            delete([paths.AnalysisID,'/MASKS_TEMPLATES/',OrigMaskName,'.HEAD']);
        end
        
        system(['3dcopy ',PLSmat.BrainMaskFile,' ',CopiedMaskFile]);
        
        % Resample mask if needed:
        ResampledMaskFile = MEGpipeline_AfniResampleAnat(CopiedMaskFile, paths.Afni4D{1}{1,1});
        
        if isempty(ResampledMaskFile)
            fprintf(ErrLog, 'ERROR: Failed to resample brainmask file: \n %s \n\n', CopiedMaskFile);
            error('ERROR: Failed to resample brainmask file for PLS.')
        end
        
        % Get brain-only indices:
        [~, BrainMask, Info, ~] = BrikLoad(ResampledMaskFile, Opt);
        
        BrainMask     = double(BrainMask);
        BrainMaskFlat = reshape(BrainMask, 1, []);       % Flatten to 1D-array (applied to flat data).
        BrainMaskFlatIndices = find(BrainMaskFlat > 0);  % Linear voxel-indices of brain-only areas.
        
        BrainMaskFlat = [];  % Free memory
        
    end
    
end  % If Brain Mask

% Why do we use linear-equivalent indices for masking? Easier & more consistent.
% - We could use the 3D-indices from the Brain/ROI Masks instead, and apply them
%   to the data (for each time) prior to flattening.
% - However, we need to flatten the data anyways, so this is just easier.



%===========================%
% ACQUIRE ROI MASK INDICES: %
%===========================%

if strcmp(PLSmat.MaskType, 'ROIMask')
	
    % Check for ROI mask:
    if ~exist(PLSmat.ROIMaskFile, 'file')
        fprintf(ErrLog, 'ERROR: Input ROI-Mask does not exist: \n %s \n\n',...
            PLSmat.ROIMaskFile);
        
        error('ERROR: Input ROI-Mask .xls could not be found.')
    end
    
    
    % Convert Mask XYZmm RAI coordinates (R-L, A-P, I-S) into AFNI IJK-indices:
    % Note: AFNI reads orientation in voxel storage-order (not +x, +y, +z).
    [ROIcoords, ROInames] = load_randy_regions(PLSmat.ROIMaskFile);
    
    for roi = 1:size(ROICoords, 1)
        [~, IJKafni] = AFNI_XYZcontinuous2Index(ROIcoords(roi,:), SampleInfo, 'RAI');
        ROIMask3DIndices(roi,:) = IJKafni + 1;  % Correct for AFNI index starting at 0.
    end
    
    
    % Get ROI voxel indices:
    ROIMaskFlatIndices = [];
    ROIBlobFlatIndices = [];
    
    for roi = 1:size(ROIMask3DIndices, 1)
        ROIBlobFlatIndices{roi} = [];
        clear Roi3DIndex RoiFlatIndex

        % Get linear-equivalent voxel-index for ROI:
        Roi3DIndex   = ROIMask3DIndices(roi,:);
        RoiFlatIndex = sub2ind(BrikDim3D, Roi3DIndex(1), Roi3DIndex(2), Roi3DIndex(3));
        ROIMaskFlatIndices = [ROIMaskFlatIndices; RoiFlatIndex];
        
        % Get indices for voxels around ROI (within ROI-blob):
        for x = PLSmat.ROIBlobSize
            for y = PLSmat.ROIBlobSize
                for z = PLSmat.ROIBlobSize
                    
                    clear BlobVoxel FlatIndex
                    BlobVoxel = Roi3DIndex + [x, y, z];  % Get 3D-index of blob-voxel.
                    
                    for dim = 1:3  % Make sure voxel is within BRIK dims.
                        BlobVoxel(dim) = max(BlobVoxel(dim), 1);
                        BlobVoxel(dim) = min(BlobVoxel(dim), BrikDim3D(dim));
                    end
                    
                    FlatIndex = sub2ind(BrikDim3D, BlobVoxel(1), BlobVoxel(2), BlobVoxel(3));
                    ROIBlobFlatIndices{roi} = [ROIBlobFlatIndices{roi}; FlatIndex];
                    
                end  % z
            end  % y
        end  % x
        
        ROIBlobFlatIndices{roi} = unique(ROIBlobFlatIndices{roi});
        
    end  % ROI
end  % If ROIMask



%=================================%
% GENERATE DATAMAT FOR PLS INPUT: %
%=================================%

% Load & flatten data, apply mask, and stack into Datamat:
% Datamat{g} for each group is stacked as Subjects within Conditions.
clear Datamat Bhvmat Seedmat BhvSeedmat

for g = 1:length(name.GroupID)
    Datamat{g} = [];
    
    for c = 1:length(name.CondID)
        for s = 1:length(name.SubjID{g})
            
            clear Brik BrikFlatXYZ BrikMasked BrikFullFlat
            [~, Brik, ~, ~] = BrikLoad(paths.Afni4D{g}{s,c}, Opt);
            Brik = double(Brik);  % Data is (x, y, z, time)

            % Flatten data for each time and apply mask:
			% Note: Cannot apply linear-equivalent of 3D-indices to 4D-data.
			% Therefore, need to apply mask separately for each time-index.
            if strcmp(PLSmat.MaskType, 'BrainMask')
                BrikFlatXYZ = reshape(Brik, [], BrikDim4D(4));      % Data is (X*Y*Z, time)
                BrikMasked  = BrikFlatXYZ(BrainMaskFlatIndices,:);  % Mask data for each time
                
            elseif strcmp(PLSmat.MaskType, 'ROIMask')
                for t = 1:BrikNumTime
                    BrikFlatXYZ = reshape(squeeze(Brik(:,:,:,t)), 1, [])  % Data is (X*Y*Z)
                    
                    % Get mean-signal of ROI-blob for each ROI:
                    for roi = 1:length(ROIBlobFlatIndices)
                        MeanBlobData      = mean(BrikFlatXYZ(ROIBlobFlatIndices{roi}));
                        BrikMasked(roi,t) = MeanBlobData;
                    end
                end
            end
            
            BrikMasked   = BrikMasked(:, [PLSmat.TimeIndices]);  % Select time-indices.
            BrikFullFlat = reshape(BrikMasked, 1, []);    % Data is 1D-array (X*Y*Z*time)
            Datamat{g}   = [Datamat{g}; BrikFullFlat];    % Add flat BRIK to group datamat.
            
        end  % Subj
    end  % Cond
	
    
    % Check datamat dimensions:
    CorrectNumRows    = (length(name.SubjID{g}) * length(name.CondID));
    
    if size(Datamat{g}, 1) ~= CorrectNumRows
        fprintf(ErrLog, ['ERROR: Datamat group has incorrect # of rows:'...
            '\n Group: %s \n Correct # rows: %s \n Current # rows: %s \n\n'],...
            name.GroupID{g}, num2str(CorrectNumRows), num2str(size(Datamat{g}, 1)));
        
        error('ERROR: Error creating Datamat. See ErrorLog.')
    end
    
end  % Group

Brik         = [];  % Free memory
BrikMasked   = [];
BrikFlatXYZ  = [];
BrikFullFlat = [];



%======================================%
% GENERATE SEEDMAT FOR BEHAVIORAL PLS: %
%======================================%
	
% Seedmat & Bhvmat arrangement:
% - Each column represents a seed or behavioral measurement.
% - Rows represent datasets (Stacked subjects within conditions within groups).
% - Subjs, Conds, and Grps should be loaded in the same order as Datamat.
% - Note: Groups are stacked prior to PLS input (unlike Datamats{g}).

if ~isempty(PLSmat.SeedFile)
    if ~exist(PLSmat.SeedFile, 'file')
        fprintf(ErrLog, 'ERROR: Seed .xls does not exist: \n %s \n\n', PLSmat.SeedFile);
        error('ERROR: Input Seed .xls file could not be found.')
    end
    
    
    % Load seed .xls data:
	xlsdata = [];
    xlsdata = importdata(PLSmat.SeedFile);
    
    if isunix || ismac  % Seed data starts at line 24.
        SeedCoords(:,1) = xlsdata.data([24:end],1);  % x
        SeedCoords(:,2) = xlsdata.data([24:end],2);  % y
        SeedCoords(:,3) = xlsdata.data([24:end],3);  % z
    elseif ispc
        SeedCoords = xlsdata.data(:,[1:3]);
    end
    
    NumSeeds    = size(SeedCoords, 1);
    SeedOrient  = xlsdata.textdata{21,2};  % orientation (voxel-order)
    SeedRegions	= xlsdata.textdata([24:end],[5:6]);
    
    if size(xlsdata.data, 2) == 3  % In case all times were specified as ranges
        xlsdata.data(:,4) = NaN;
    end
    
    
    % Get time-indices for each seed-coordinate:
    TempDataCoords(:,1) = xlsdata.data([24:end],4);      % time
    TempTextCoords(:,1) = xlsdata.textdata([24:end],4);  % in case where time-range specified
    
    for seed = 1:NumSeeds
        if ~isnan(TempDataCoords(seed,1))
            SeedTimes{seed} = TempDataCoords(seed,1);
        else
            SeedTimes{seed} = str2num(TempTextCoords{seed});
            
            if isempty(SeedTimes{seed})
                fprintf(ErrLog, 'ERROR: Time-indices for seeds were not entered properly. \n\n');
                error('ERROR: Time-indices for seeds were not entered properly.')
            end
        end
    end
    
    SeedTimes = SeedTimes';
    clear TempDataCoords TempTextCoords
    
    
	% Check seed .xls data:
    if ~strcmp(xlsdata.textdata{21,1}, 'Orientation:')
        fprintf(ErrLog, 'ERROR: File is not a proper seed .xls file. \n\n');
        error('ERROR: File loaded is not a proper seed .xls file.')
    end
    
    if ~isempty(find(isnan(SeedCoords)))
        fprintf(ErrLog, 'ERROR: SeedCoords contain NaNs (Missing data cells?). \n\n');
        error('ERROR: SeedCoord data contains NaNs. Check .xls for missing cells.')
    end
    
    if ~isempty(find(cellfun(@isempty, SeedRegions)))
        fprintf(ErrLog, 'ERROR: Some seeds are missing region names. \n\n');
        error('ERROR: Seeds are missing region names. Check .xls file.')
    end
    
    if isempty(SeedOrient)
        fprintf(ErrLog, 'ERROR: SeedCoord orient not specified in .xls file. \n\n');
        error('ERROR: SeedCoord orientation not specified in .xls file.')
    end
    
    AllTimeIndices = 1:size(time.Windows, 1);
    for seed = 1:NumSeeds
        if min(ismember(SeedTimes{seed}, AllTimeIndices)) == 0
            fprintf(ErrLog, 'ERROR: Seed time-indices exceed data time-indices. \n\n');
            error('ERROR: Seed time-indices exceed data time-indices.')
        end
    end
    
    
    % Generate name for each seed:
    for seed = 1:NumSeeds
        SeedTimeIndex = SeedTimes{seed};
        
        if length(SeedTimeIndex) == 1
            SeedNames{seed} = [SeedRegions{seed,1},'_',...
                SeedRegions{seed,2},'_',time.str.Windows{SeedTimeIndex,3}];
        else
            SeedNames{seed} = [SeedRegions{seed,1},'_',...
                SeedRegions{seed,2},'_Avg',num2str(length(SeedTimeIndex)),'TimeIndices'];
        end
        
        SeedNames{seed} = deblank(SeedNames{seed});
        SeedNames{seed} = regexprep(SeedNames{seed}, '\s+', '_');
    end
    
    
    % Convert seed coords from XYZmm to AFNI 4D IJK-index:
    % Note: AFNI reads orientation in voxel storage-order (not +x, +y, +z).
    if strcmp(SeedOrient, 'IJK') || strcmp(SeedOrient, 'IJKafni')
        Seeds3DIndices = SeedCoords;
        
    else
        for seed = 1:NumSeeds
            clear SeedXYZ IJKafni
            
            SeedXYZ      = SeedCoords(seed,:);  % XYZmm of current seed
            [~, IJKafni] = AFNI_XYZcontinuous2Index(SeedXYZ, SampleInfo, SeedOrient);
            Seeds3DIndices(seed,:) = [IJKafni+1];  % (I, J, K of current seed)
        end
    end
    
    
    % Get voxel indices for each seed:
    SeedFlatIndices     = [];
    SeedBlobFlatIndices = [];
    
    for seed = 1:NumSeeds
        SeedFlatIndices{seed}     = [];
        SeedBlobFlatIndices{seed} = [];
        clear Seed3DIndex
		
        % Get linear-equivalent voxel-index for seed:
        Seed3DIndex = Seeds3DIndices(seed,:);
        
        for t = SeedTimes{seed}  % see below for description
            clear SeedFlatIndex
            
            SeedFlatIndex = sub2ind...
                (BrikDim4D, Seed3DIndex(1), Seed3DIndex(2), Seed3DIndex(3), t);
            
            SeedFlatIndices{seed} = [SeedFlatIndices{seed}; SeedFlatIndex];
            
            % Get linear-equivalent indices for voxels around seed (within seed-blob):
            % Later on, data from all voxels within the "seed blob" will be averaged.
            for x = PLSmat.SeedBlobSize
                for y = PLSmat.SeedBlobSize
                    for z = PLSmat.SeedBlobSize
                        
                        clear BlobVoxel FlatIndex
                        BlobVoxel = [Seed3DIndex,t] + [x, y, z, 0];  % Get 4D-index of blob-voxel.
                        
                        for dim = 1:4  % Make sure voxel is within BRIK dims.
                            BlobVoxel(dim) = max(BlobVoxel(dim), 1);
                            BlobVoxel(dim) = min(BlobVoxel(dim), BrikDim4D(dim));
                        end
                        
                        FlatIndex = sub2ind(BrikDim4D,...
                            BlobVoxel(1), BlobVoxel(2), BlobVoxel(3), BlobVoxel(4));
                        
                        SeedBlobFlatIndices{seed} = [SeedBlobFlatIndices{seed}; FlatIndex];
                        
                    end  % z
                end  % y
            end  % x
        end  % t
        
        SeedBlobFlatIndices{seed} = unique(SeedBlobFlatIndices{seed});
        
        % Later on, for each seed, all indices within the "seed blob" will be averaged.
        % - In the case of users specifying a time-range or multiple time-indices for
        %   a given seen, we add voxel-data at all specified time-indices to the seed blob.
        % - Since all data from all voxel indices within a "seed blob" are averaged later
        %   on, this will mean data from each time-index is also averaged as desired.
        
    end  % Seed
	
    
	% Extract data for each seed and load into {Group}(Subj,Cond,Bhv):
	GroupSeedData = [];
    
	for g = 1:length(name.GroupID)
		for s = 1:length(name.SubjID{g})
			for c = 1:length(name.CondID)
                
				[~, Brik, Info, ~] = BrikLoad(paths.Afni4D{g}{s,c}, Opt);
				Brik = double(Brik);

				% Note: Do not need to flatten data prior to applying seed-indices.
				% Applying linear-equivalent of 4D-indices to 4D-data.
				for seed = 1:NumSeeds
					GroupSeedData{g}(s,c,seed) = mean(Brik(SeedBlobFlatIndices{seed}));
				end
				
			end  % Cond
		end  % Subj
	end  % Group
	
    
	% Reshape to stack Subjects within Conditions: {Group}(Subj*Cond,Seed)
	% Then stack Groups to get Seedmat (Subj*Cond*Group, Seed).
	Seedmat = [];
	for g = 1:length(name.GroupID)
		GroupSeedData{g} = reshape(GroupSeedData{g}, [], NumSeeds);
		Seedmat          = [Seedmat; GroupSeedData{g}];
    end
	
    
	% Check dims on Seedmat:
    CorrectNumRows = sum(cellfun(@numel, name.SubjID)) * length(name.CondID);
    
    if size(Seedmat, 1) ~= CorrectNumRows
        fprintf(ErrLog, ['ERROR: Seedmat has incorrect # of rows:'...
            '\n Correct # Rows: %s \n Current # rows: %s \n\n'],...
            num2str(CorrectNumRows), num2str(size(Seedmat, 1)));
        
        error('ERROR: Error creating Seedmat. See ErrorLog.')
    end
    
    Brik          = [];  % Free memory
    
end  % If SeedFile



%===========================================%
% GENERATE BEHAVIOR MAT FOR BEHAVIORAL PLS: %
%===========================================%

% Seedmat & Bhvmat arrangement:
% - Each column represents a seed or behavioral measurement.
% - Rows represent datasets (Stacked subjects within conditions within groups).
% - Subjs, Conds, and Grps should be loaded in the same order as Datamat.
% - Note: Groups are stacked prior to PLS input (unlike Datamats{g}).

if ~isempty(PLSmat.BhvFile)
    if ~exist(PLSmat.BhvFile, 'file')
        fprintf(ErrLog, 'ERROR: Behavior .xls does not exist: \n %s \n\n', PLSmat.BhvFile);
        error('ERROR: Input behavior .xls file could not be found.')
    end
    
    
	% Load behavior file:
	xlsdata = [];
    xlsdata = importdata(PLSmat.BhvFile);

    BhvSheets = fieldnames(xlsdata.textdata);
    BhvSheets(strcmp(BhvSheets, 'Instructions')) = [];
    NumBhvs   = length(BhvSheets);
    
    TotalNumRows = sum(cellfun(@numel, name.SubjID));  % #Rows = Total #Subj
    TotalNumCols = length(name.CondID);
    
    
    % Check if .xls file is proper format:
    if ~strcmp(xlsdata.textdata.(BhvSheets{1}){1,1}, 'Behavior:')
        fprintf(ErrLog, 'ERROR: File is not a proper behavior .xls file. \n\n');
        error('ERROR: File loaded is not a proper behavior .xls file.')
    end
	
    
	% Get behavior sheets & check for errors:
	for b = 1:NumBhvs
        BhvData{b}  = xlsdata.data.(BhvSheets{b});
        BhvNames{b} = xlsdata.textdata.(BhvSheets{b}){1,2};
        BhvNames{b} = deblank(BhvNames{b});
        BhvNames{b} = regexprep(BhvNames{b}, '\s+', '_');
           
		if size(BhvData{b}, 1) ~= TotalNumRows  % Rows = Group+Subj in bhv .xls
			fprintf(ErrLog, ['ERROR: Behavior .xls file has incorrect # of rows:'...
				'\n Correct # Rows = Total # Subj: %s \n Current # Rows: %s \n\n'],...
				num2str(TotalNumRows), num2str(size(BhvData{b}, 1)));
            
			error('ERROR: Behavior .xls file has incorrect # rows. See ErrorLog.')
        end
        
		if size(BhvData{b}, 2) ~= TotalNumCols  % Columns = Conds in bhv .xls
			fprintf(ErrLog, ['ERROR: Behavior .xls file has incorrect # of columns:'...
				'\n Correct # Cols = # Conds: %s \n Current # Cols: %s \n\n'],...
				num2str(TotalNumCols), num2str(size(BhvData{b}, 2)));
            
			error('ERROR: Behavior .xls file has incorrect # columns. See ErrorLog.')
        end
        
		if ~isempty(find(isnan(BhvData{b})))
			fprintf(ErrLog, 'ERROR: Behavior data contains NaNs (Missing cells?) \n\n');
			error('ERROR: Behavior data contains NaNs. Check .xls for missing cells.');
        end
        
		if isempty(BhvNames{b})
			fprintf(ErrLog, 'ERROR: Behavior name is missing. \n\n');
			error('ERROR: No name specified for behavior.');
        end
        
	end  % Bhv

    
	% Load behavior data into {Group}(Subj,Cond,Bhv) for easy reshaping:
	% For each group, get data from corresponding subject rows:
	GroupBhvData = [];
	for b = 1:NumBhvs
		StartRow = 1;  % Reset for new bhv sheet
        
		for g = 1:length(name.GroupID)
			EndRow                 = StartRow + length(name.SubjID{g}) - 1;
			GroupBhvData{g}(:,:,b) = BhvData{b}([StartRow:EndRow],[1:end]);
			StartRow               = StartRow + length(name.SubjID{g});
		end
    end
	
    
	% Reshape to stack Subjects within Conditions: {Group}(Subj*Cond,Bhv)
	% Then stack Groups to get BehaviorMat (Subj*Cond*Group, Bhv).
	Bhvmat = [];
	for g = 1:length(name.GroupID)
		GroupBhvData{g} = reshape(GroupBhvData{g}, [], NumBhvs);
		Bhvmat          = [Bhvmat; GroupBhvData{g}];
    end
	
    
    % Check dims on Bhvmat:
	CorrectNumRows = sum(cellfun(@numel, name.SubjID)) * length(name.CondID);
    
    if size(Bhvmat, 1) ~= CorrectNumRows
        fprintf(ErrLog, ['ERROR: Bhvmat has incorrect # of rows:'...
            '\n Correct # Rows: %s \n Current # rows: %s \n\n'],...
            num2str(CorrectNumRows), num2str(size(Bhvmat, 1)));
        
        error('ERROR: Error creating Bhvmat. See ErrorLog.')
    end
	
end  % If BhvFile
    


%=========================================%
% IF BOTH LOADED, MERGE BHVMAT & SEEDMAT: %
% Note: Seeds are loaded after Behaviors. %
%=========================================%

if ~isempty(PLSmat.SeedFile) && ~isempty(PLSmat.BhvFile)
	BhvSeedmat   = [Bhvmat, Seedmat];
	BhvSeedNames = [BhvNames, SeedNames];
	
    if size(BhvSeedmat, 2) ~= (NumBhvs + NumSeeds)
        fprintf(ErrLog, ['ERROR: BhvSeedmat has incorrect # of columns.:'...
            '\n Correct # Cols = NumBhvs + NumSeeds = %s \n Current # Cols: %s \n\n'],...
            num2str(NumBhvs + NumSeeds), num2str(size(BhvSeedmat, 2)));
        
        error('ERROR: Error merging behaviors and seeds. See ErrorLog.')
    end
    
end



%============================%
% SAVE VARIABLES TO PLS MAT: %
%============================%

% Get datamat & brik info:
PLSmat.Datamat        = Datamat;
PLSmat.BrikDim3D      = BrikDim3D;
PLSmat.BrikDim4D      = BrikDim4D;
PLSmat.BrikNumTime    = BrikDim4D(4);
PLSmat.BrikOrient	  = Info.ORIENT_SPECIFIC;
PLSmat.BrikSampleInfo = SampleInfo;


% Get mask info:
switch PLSmat.MaskType
	case 'BrainMask'
        PLSmat.BrainMask            = BrainMask;
		PLSmat.BrainMaskFlatIndices = BrainMaskFlatIndices;
        
	case 'ROIMask'
		PLSmat.ROIMaskCoords.RAImm = ROIcoords;
		PLSmat.ROIMaskCoords.Names = ROInames;
		PLSmat.ROIMask3DIndices    = ROIMask3DIndices;
		PLSmat.ROIMaskFlatIndices  = ROIMaskFlatIndices;
		PLSmat.ROIBlobsFlatIndices = ROIBlobsFlatIndices;
end


% Get bhvmat info:
if exist('BhvSeedmat', 'var') && ~isempty(PLSmat.BhvFile) && ~isempty(PLSmat.SeedFile)
	PLSmat.BhvType    = 'BhvSeed';
	PLSmat.Bhvmat     = BhvSeedmat;
	PLSmat.BhvNames   = BhvSeedNames;
	PLSmat.SeedCoords = SeedCoords;
	PLSmat.SeedOrient = SeedOrient;
    
elseif exist('Bhvmat', 'var') && ~isempty(PLSmat.BhvFile)
	PLSmat.BhvType  = 'Bhv';
	PLSmat.Bhvmat   = Bhvmat;
	PLSmat.BhvNames = BhvNames;
    
elseif exist('Seedmat', 'var') && ~isempty(PLSmat.SeedFile)
	PLSmat.BhvType    = 'Seed';
	PLSmat.Bhvmat     = Seedmat;
	PLSmat.BhvNames   = SeedNames;
	PLSmat.SeedCoords = SeedCoords;
    PLSmat.SeedTimes  = SeedTimes;
	PLSmat.SeedOrient = SeedOrient;
	%** save SeedsFlatIndices and SeedBlobsFlatIndices?
end



%=================%

if exist([pwd,'/ErrorLog_MakePLSmat.txt'], 'file')
    LogCheck = dir('ErrorLog_MakePLSmat.txt');
    if LogCheck.bytes ~= 0  % File is not empty
        open('ErrorLog_MakePLSmat.txt');
    else
        delete('ErrorLog_MakePLSmat.txt');
    end
end

fclose(ErrLog);
diary off